Load data

source("~/Dropbox/research/atacseq/method/gcqn_validated.R")
library(Rsubread)
library(GenomicAlignments)
## Loading required package: BiocGenerics
## Loading required package: parallel
## 
## Attaching package: 'BiocGenerics'
## The following objects are masked from 'package:parallel':
## 
##     clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
##     clusterExport, clusterMap, parApply, parCapply, parLapply,
##     parLapplyLB, parRapply, parSapply, parSapplyLB
## The following objects are masked from 'package:stats':
## 
##     IQR, mad, sd, var, xtabs
## The following objects are masked from 'package:base':
## 
##     anyDuplicated, append, as.data.frame, basename, cbind, colnames,
##     dirname, do.call, duplicated, eval, evalq, Filter, Find, get, grep,
##     grepl, intersect, is.unsorted, lapply, Map, mapply, match, mget,
##     order, paste, pmax, pmax.int, pmin, pmin.int, Position, rank,
##     rbind, Reduce, rownames, sapply, setdiff, sort, table, tapply,
##     union, unique, unsplit, which, which.max, which.min
## Loading required package: S4Vectors
## Loading required package: stats4
## 
## Attaching package: 'S4Vectors'
## The following object is masked from 'package:base':
## 
##     expand.grid
## Loading required package: IRanges
## Loading required package: GenomeInfoDb
## Loading required package: GenomicRanges
## Loading required package: SummarizedExperiment
## Loading required package: Biobase
## Welcome to Bioconductor
## 
##     Vignettes contain introductory material; view with
##     'browseVignettes()'. To cite Bioconductor, see
##     'citation("Biobase")', and for packages 'citation("pkgname")'.
## Loading required package: DelayedArray
## Loading required package: matrixStats
## 
## Attaching package: 'matrixStats'
## The following objects are masked from 'package:Biobase':
## 
##     anyMissing, rowMedians
## Loading required package: BiocParallel
## 
## Attaching package: 'DelayedArray'
## The following objects are masked from 'package:matrixStats':
## 
##     colMaxs, colMins, colRanges, rowMaxs, rowMins, rowRanges
## The following objects are masked from 'package:base':
## 
##     aperm, apply, rowsum
## Loading required package: Biostrings
## Loading required package: XVector
## 
## Attaching package: 'Biostrings'
## The following object is masked from 'package:base':
## 
##     strsplit
## Loading required package: Rsamtools
library(scales)
library(tidyverse)
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.0 ──
## ✔ ggplot2 3.3.2     ✔ purrr   0.3.3
## ✔ tibble  2.1.3     ✔ dplyr   0.8.4
## ✔ tidyr   1.0.2     ✔ stringr 1.4.0
## ✔ readr   1.3.1     ✔ forcats 0.4.0
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ readr::col_factor() masks scales::col_factor()
## ✖ dplyr::collapse()   masks Biostrings::collapse(), IRanges::collapse()
## ✖ dplyr::combine()    masks Biobase::combine(), BiocGenerics::combine()
## ✖ purrr::compact()    masks XVector::compact()
## ✖ dplyr::count()      masks matrixStats::count()
## ✖ dplyr::desc()       masks IRanges::desc()
## ✖ purrr::discard()    masks scales::discard()
## ✖ tidyr::expand()     masks S4Vectors::expand()
## ✖ dplyr::filter()     masks stats::filter()
## ✖ dplyr::first()      masks GenomicAlignments::first(), S4Vectors::first()
## ✖ dplyr::lag()        masks stats::lag()
## ✖ dplyr::last()       masks GenomicAlignments::last()
## ✖ ggplot2::Position() masks BiocGenerics::Position(), base::Position()
## ✖ purrr::reduce()     masks GenomicRanges::reduce(), IRanges::reduce()
## ✖ dplyr::rename()     masks S4Vectors::rename()
## ✖ purrr::simplify()   masks DelayedArray::simplify()
## ✖ dplyr::slice()      masks XVector::slice(), IRanges::slice()
library(ggplot2)
library(qsmooth)
library(edgeR)
## Loading required package: limma
## 
## Attaching package: 'limma'
## The following object is masked from 'package:BiocGenerics':
## 
##     plotMA
library(cqn)
## Loading required package: mclust
## Package 'mclust' version 5.4.5
## Type 'citation("mclust")' for citing this R package in publications.
## 
## Attaching package: 'mclust'
## The following object is masked from 'package:purrr':
## 
##     map
## Loading required package: nor1mix
## Loading required package: preprocessCore
## Loading required package: splines
## Loading required package: quantreg
## Loading required package: SparseM
## 
## Attaching package: 'SparseM'
## The following object is masked from 'package:base':
## 
##     backsolve
gcqn_qsmooth <- function(counts, gcGroups, bio){
  gcBinNormCounts <- matrix(NA, nrow=nrow(counts), ncol=ncol(counts), dimnames=list(rownames(counts),colnames(counts)))
  for(ii in 1:nlevels(gcGroups)){
    id <- which(gcGroups==levels(gcGroups)[ii])
    countBin <- counts[id,]
    qs <- qsmooth(countBin, group_factor=bio)
    normCountBin <- qs@qsmoothData
    normCountBin <- round(normCountBin)
    normCountBin[normCountBin<0] <- 0
    gcBinNormCounts[id,] <- normCountBin
  }
  return(gcBinNormCounts)
}
# note that genome patch is likely not identical to the one used by Tiffany for mapping.
# download.file("ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCA/000/001/635/GCA_000001635.8_GRCm38.p6/GCA_000001635.8_GRCm38.p6_genomic.fna.gz", destfile="/home/compomics/koen/genomes/mm10/p6/GCA_000001635.8_GRCm38.p6_genomic.fna.gz")

makeGR_broadPeak <- function(file, genome){
  ## this is for broadPeak files, e.g. from MACS2.

  # import peaks and make GenomicRanges
  peaks <- read.delim(file, header=FALSE)
  peakID <- unlist(lapply(strsplit(as.character(peaks[,4]),split="_"),function(x) paste(x[4:5], collapse="_")))
  peaksGr <- GRanges(seqnames=gsub(x=as.character(peaks[,1]), pattern="chr", replacement=""), ranges=IRanges(start=peaks[,2], end=peaks[,3]), strand=NA, mcols=DataFrame(peakID=peakID))
  names(peaksGr) <- mcols(peaksGr)[,1]

  # get GC content of peaks
  ff <- FaFile(genome)
  peakSeqs <- getSeq(x=ff, peaksGr)
  gcContentPeaks <- letterFrequency(peakSeqs, "GC",as.prob=TRUE)[,1]
  mcols(peaksGr) <- cbind(mcols(peaksGr), gc=gcContentPeaks)
  return(peaksGr)
}

# MACS paired-end
grMacsPE <- makeGR_broadPeak(file="~/Dropbox/compomicsVM/protocolDAta/pooledPeaksMACS2_allProtocolSamples_cutoff10_peaks.broadPeak", genome="~/data/genomes/mouse/Mus_musculus.GRCm38.dna.primary_assembly.fa")
## Warning in .normarg_strand(strand, seqnames): missing values in 'strand'
## converted to "*"
fcMacsPE <- readRDS("~/Dropbox/compomicsVM/protocolDAta/featureCounts_MACS2BAMPE.rds")


# get metadata
library(openxlsx)
metadata <- read.xlsx("~/Dropbox/research/atacseq/bulk/analysis/protocolData/20180611_atacseq_metadata.xlsx")
index <-  metadata[,"Index.#"]
batch <- as.factor(ifelse(metadata$Transposition.date==43228,"A","B"))
cellNumber <- metadata$Cell.number
atacMethod <- as.factor(metadata$ATACseq.method)
pcrCycle <- metadata$PCR.cycles
cellType <- metadata$Cell.type

knitr::kable(metadata[,2:5])
Animal.&.treatment Cell.type Cell.number ATACseq.method
2 male CD-1, DOB 4/12/18, 3.7w old sorted ICAM1- cells from dissected OE 5000 original
2 male CD-1, DOB 4/12/18, 3.7w old sorted ICAM1- cells from dissected OE 10000 original
2 male CD-1, DOB 4/12/18, 3.7w old sorted ICAM1- cells from dissected OE 50000 original
2 male CD-1, DOB 4/12/18, 3.7w old sorted ICAM1+ cells (mix of PE+ or FITC+) from dissected OE 10355 original
2 male CD-1, DOB 4/12/18, 3.7w old sorted ICAM1+ cells (mix of PE+ or FITC+) from dissected OE 10355 Omni
1 male CD-1, DOB 4/12/18, 4.7w old sorted ICAM1- cells from dissected OE 10000 Omni
1 male CD-1, DOB 4/12/18, 4.7w old sorted ICAM1- cells from dissected OE 10000 Omni
1 male CD-1, DOB 4/12/18, 4.7w old sorted ICAM1- cells from dissected OE 10000 Omni
1 male CD-1, DOB 4/12/18, 4.7w old sorted ICAM1+ cells from dissected OE 6053 original
2 male CD-1, DOB 4/12/18, 4.7w old, 24 HPI of methimazole sorted ICAM1+ cells from dissected OE 15154 original

EDA

plotWidthHex <- function(gr, counts){
  counts2 <- counts
  colnames(counts2) <- paste0("Batch_",batch,",n=",cellNumber,",method_",atacMethod,",PCRcycs_",pcrCycle,",id",index)
  df <- as_tibble(cbind(counts2,width=width(gr)))
  df <- gather(df, sample, value, -width)
  ggplot(data=df, aes(x=log(width+1), y=log(value+1)) ) + ylab("log(count + 1)") + xlab("log(width + 1)") + geom_hex(bins = 50) + theme_bw() + facet_wrap(~sample, nrow=2) #+ annotate("text",label=cellType, x=0.5, y=12, size=3)
}
plotWidthHex(gr=grMacsPE, fcMacsPE$counts) + xlim(3,10) + ggtitle("Peak caller: MACS PE")

plotGCHex <- function(gr, counts){
  counts2 <- counts
  colnames(counts2) <- paste0("Batch_",batch,",n=",cellNumber,",method_",atacMethod,",PCRcycs_",pcrCycle,",id",index)
  df <- as_tibble(cbind(counts2,gc=mcols(gr)$gc))
  df <- gather(df, sample, value, -gc)
  ggplot(data=df, aes(x=gc, y=log(value+1)) ) + ylab("log(count + 1)") + xlab("GC content") + geom_hex(bins = 50) + theme_bw() + facet_wrap(~sample, nrow=2) #+ annotate("text",label=cellType, x=0.5, y=12, size=3)
}
plotGCHex(gr=grMacsPE, fcMacsPE$counts) + ggtitle("Peak caller: MACS PE")

evaluateSimulation <- function(simCounts, simGC, simWidth, design, deId, grSim, nSamples, grpSim){
  testEdgeR <- function(counts, design, tmm=FALSE, offset=NULL){
    d <- DGEList(counts)
    if(tmm) d <- calcNormFactors(d)
    if(!is.null(offset)) d$offset <- offset
    d <- estimateDisp(d, design)
    fit <- glmFit(d, design)
    lrt <- glmLRT(fit, coef=2)
    lrt$table$padj <- p.adjust(lrt$table$PValue, "fdr")
    return(lrt)
  }
  
  library(DESeq2)
  testDESeq2 <- function(counts, covarDf){
    dds <- DESeqDataSetFromMatrix(counts,
                                colData=covarDf,
                                design=as.formula("~ group"))
    dds <- DESeq(dds)
    res <- results(dds)
    return(res)
  }
  
  ## normalize
  ## TMM
  d <- DGEList(simCounts)
  d <- calcNormFactors(d)
  cpmTMM <- edgeR::cpm(d, normalized.lib.sizes=TRUE, log=TRUE, prior.count=0.25)
  
  ## FQN
  fqCounts <- FQnorm(simCounts)
  cpmFQ <- edgeR::cpm(fqCounts, normalized.lib.sizes=FALSE, log=TRUE, prior.count=0.25)
  
  ## qsmooth
  qsmoothObj <- qsmooth(object = simCounts, group_factor = grpSim)
  qsmoothCounts <- qsmoothData(qsmoothObj)
  
  ## cqn
  cqnObj <- cqn(simCounts, x=simGC, lengths=simWidth, sizeFactors=colSums(simCounts))
  cqnplot(cqnObj, xlab="GC content", lwd=2)
  cqnCounts <- 2^(cqnObj$y + cqnObj$offset)
  cpmCqn <- edgeR::cpm(cqnCounts, normalized.lib.sizes=FALSE, log=TRUE, prior.count=0.25)
  
  ## GCQN:
  gcGroups <- Hmisc::cut2(simGC, g=50)
  countsGCQN <- gcqn(simCounts, gcGroups, summary="median", round=TRUE)
  cpmGCQN <- edgeR::cpm(countsGCQN, normalized.lib.sizes=FALSE, log=TRUE, prior.count=0.25)
  
  ## smooth GCQN
  normCountsGcqnSmooth <- gcqn_qsmooth(simCounts, gcGroups, bio=grpSim)
  
  ## EDASeq
  library(EDASeq)
  wit <- withinLaneNormalization(x=as.matrix(simCounts), y=simGC, which="full")
  countsEDASeq <- betweenLaneNormalization(wit, "full")
  cpmEDASeq <- edgeR::cpm(countsEDASeq, normalized.lib.sizes=FALSE, log=TRUE, prior.count=0.25)
  
  ## RUVg
  library(RUVSeq)
  hk <- readRDS("~/data/genomes/hkListMouseGenomicRanges.rds")
  qh <- queryHits(findOverlaps(grSim[-deId], hk, type="within"))
  #qh <- queryHits(findOverlaps(grSim, hk, type="within"))
  negcon <- rownames(simCounts)[qh]
  ruvRes <- RUVg(as.matrix(simCounts), negcon, k=2)
  ruvgCounts <- ruvRes$normalizedCounts
  cpmRUV <- edgeR::cpm(ruvgCounts, normalized.lib.sizes=FALSE, log=TRUE, prior.count=0.25)
  
  ## RUVs
  # scid <- matrix(c(which(design[,2] == 0),
  #          which(design[,2] == 1)),
  #        nrow=2, ncol=nSamples/2, byrow=TRUE)
  # ruvsRes <- RUVs(as.matrix(simCounts), 
  #                 cIdx = 1:nrow(simCounts),
  #                 scIdx = scid,
  #                 k=2)
  # ruvsCounts <- ruvsRes$normalizedCounts

  
  
  # TMM
  resTMM <- testEdgeR(simCounts, design, tmm=TRUE)
  # DESeq2 MOR
  resDESeq2 <- testDESeq2(simCounts, covarDf=data.frame(group=grpSim))
  # Cqn
  resCqn <- testEdgeR(simCounts, design, tmm=FALSE, offset=cqnObj$glm.offset)
  #FQ
  resFQ <- testEdgeR(fqCounts, design, tmm=FALSE)
  # qsmooth
  resQsmooth <- testEdgeR(qsmoothCounts, design, tmm=FALSE)
  # EDASEQ
  resEDASeq <- testEdgeR(countsEDASeq, design, tmm=FALSE)
  # GCQN
  resGCQN <- testEdgeR(countsGCQN, design, tmm=FALSE)
  # smooth GCQN
  resGCQNSmooth <- testEdgeR(normCountsGcqnSmooth, design, tmm=FALSE)
  # RUVg
  resRUVg <- testEdgeR(ruvgCounts, design, tmm=FALSE)
  # RUVs
  #resRUVs <- testEdgeR(ruvsCounts, design, tmm=FALSE)

  
  truth <- rep(0, nrow(simCounts))
  truth[deId] <- 1
  truthDf <- data.frame(truth=truth)
  rownames(truthDf) <- rownames(simCounts)
  
  library(iCOBRA)
  cbd1 <- COBRAData(pval = data.frame(TMM=resTMM$table$PValue,
                                     DESeq2=resDESeq2$pvalue,
                                     qsmooth=resQsmooth$table$PValue,
                                      cqn=resCqn$table$PValue,
                                      FQ=resFQ$table$PValue,
                                      edaseq=resEDASeq$table$PValue,
                                      gcqn=resGCQN$table$PValue,
                                      gcqn_smooth=resGCQNSmooth$table$PValue,
                                      RUVg=resRUVg$table$PValue,
                                      #RUVs=resRUVs$table$PValue,
                                      row.names=rownames(simCounts)),
            truth = truthDf)
  # cbd <- calculate_adjp(cbd1)
  # cbd <- calculate_performance(cbd, binary_truth = "truth")
  # cbp <- prepare_data_for_plot(cbd)
  # p <- plot_fdrtprcurve(cbp, pointsize=2, yaxisrange=c(0,1),
  #                       title=paste0(nSamples," samples"))
  # print(p)
  #return(list(cbd=cbd1, p=p))
  return(cbd1)
}


simulateAndEvaluate <- function(seed, nTotal, inputCounts, grp, nDE, meanlog=0.8, sdlog=0.1, simGC, simWidth, grSim){
  
  set.seed(seed)
  n1 <- 1:(nTotal/2)
  n2 <- (nTotal/2+1):nTotal
  fractions <- sweep(inputCounts, 2, colSums(inputCounts), "/")
  inputCountssim <- matrix(NA, nrow=nrow(inputCounts), ncol=nTotal,
                     dimnames=list(rownames(inputCounts), paste0("sample", 1:nTotal)))
  remainingSamplesA <- which(grp == "A")
  remainingSamplesB <- which(grp == "B")
  
  deId <- sample(x=1:nrow(inputCounts), size=nDE)
  delta <- rep(0, nrow(inputCounts))
  sign <- rep(c(1,-1), each=nDE/2)
  if(length(sign) != nDE) sign <- c(sign,1)
  # delta[deId] <- sign*log(rlnorm(n=1e4, meanlog=3, sdlog=1.2))
  delta[deId] <- sign*log(rlnorm(n=nDE, meanlog=meanlog, sdlog=sdlog))
  hist(delta[deId],  breaks=40)
  
  # group 1
  for(ss in n1){
    sampleId <- sample(x=remainingSamplesA, size=1)
    remainingSamplesA <- remainingSamplesA[!remainingSamplesA %in% sampleId]
    curLS <- runif(n=1, min(colSums(inputCounts)), max(colSums(inputCounts)))
    curFracs <- fractions[,sampleId]
    curFracs[delta < 0] <- (curFracs[delta < 0]) * exp(abs(delta[delta < 0]))
    simY <- rmultinom(n=1, size=curLS, prob=curFracs)
    inputCountssim[,ss] <- simY
  }
  
  # group 2
  for(ss in n2){
    sampleId <- sample(x=remainingSamplesB, size=1)
    remainingSamplesB <- remainingSamplesB[!remainingSamplesB %in% sampleId]
    curLS <- runif(n=1, min(colSums(inputCounts)), max(colSums(inputCounts)))
    curFracs <- fractions[,sampleId]
    curFracs[delta > 0] <- (curFracs[delta > 0]) * exp(abs(delta[delta > 0]))
    simY <- rmultinom(n=1, size=curLS, prob=curFracs)
    inputCountssim[,ss] <- simY
  }
  
  plot(x=delta, y=log(rowMeans(inputCountssim[,n2]) / rowMeans(inputCountssim[,n1])))
  abline(0,1, col="red")
  abline(h=0, col="blue", lty=2)
  
  grpSim <- factor(rep(c("A", "B"), each=nTotal/2))
  evaluateSimulation(simCounts = inputCountssim, 
                     simGC = simGC, 
                     simWidth = simWidth, 
                     design = model.matrix(~grpSim), 
                     deId = deId, 
                     grSim = grSim,
                     nSamples = nTotal,
                     grpSim = grpSim)
}
inputCounts <- fcMacsPE$counts
ct <- metadata$Cell.type
# grp <- vector(length=length(ct))
# grp[grep(x=ct, pattern="ICAM1-", fixed=TRUE)] <- "A"
# grp[grep(x=ct, pattern="ICAM1+", fixed=TRUE)] <- "B"
simGC <- mcols(grMacsPE)$gc
simWidth <- width(grMacsPE)

nDEParams <- round(nrow(inputCounts)/c(10,4))
nTotalParams <- c(8, 10)
for(nde in 1:length(nDEParams)){
  for(ntot in 1:length(nTotalParams)){
      resList <- list()
    for(iter in 1:14){
      grp <- sample(rep(c("A", "B"), each=5))
      resList[[iter]] <- simulateAndEvaluate(seed=iter, 
                                             nTotal=nTotalParams[ntot], 
                                             inputCounts=inputCounts, 
                                             grp=grp, 
                                             nDE=nDEParams[nde], 
                                             meanlog=0.8, 
                                             sdlog=0.1, 
                                             simGC=simGC, 
                                             simWidth=simWidth, 
                                             grSim=grMacsPE)
    }
    saveRDS(resList, file=paste0("simRes_n",nTotalParams[ntot],"_nDE",nde,".rds"))
  }
}

## Warning: The use of 'sig2' is deprecated; do specify 'sigma' (= sqrt(sig2))
## instead

## Loading required package: ShortRead
## 
## Attaching package: 'ShortRead'
## The following object is masked from 'package:dplyr':
## 
##     id
## The following object is masked from 'package:purrr':
## 
##     compose
## The following object is masked from 'package:tibble':
## 
##     view
## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
## 
## Attaching package: 'iCOBRA'
## The following object is masked from 'package:Biostrings':
## 
##     score
## The following objects are masked from 'package:GenomicRanges':
## 
##     score, score<-
## The following objects are masked from 'package:BiocGenerics':
## 
##     score, score<-

## Warning: The use of 'sig2' is deprecated; do specify 'sigma' (= sqrt(sig2))
## instead

## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
## Warning: Zero sample variances detected, have been offset away from zero

## Warning: The use of 'sig2' is deprecated; do specify 'sigma' (= sqrt(sig2))
## instead

## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
## Warning: Zero sample variances detected, have been offset away from zero

## Warning: The use of 'sig2' is deprecated; do specify 'sigma' (= sqrt(sig2))
## instead

## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing

## Warning: The use of 'sig2' is deprecated; do specify 'sigma' (= sqrt(sig2))
## instead

## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing

## Warning: The use of 'sig2' is deprecated; do specify 'sigma' (= sqrt(sig2))
## instead

## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing

## Warning: The use of 'sig2' is deprecated; do specify 'sigma' (= sqrt(sig2))
## instead

## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing

## Warning: The use of 'sig2' is deprecated; do specify 'sigma' (= sqrt(sig2))
## instead

## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
## Warning: Zero sample variances detected, have been offset away from zero

## Warning: The use of 'sig2' is deprecated; do specify 'sigma' (= sqrt(sig2))
## instead

## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing

## Warning: The use of 'sig2' is deprecated; do specify 'sigma' (= sqrt(sig2))
## instead

## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing

## Warning: The use of 'sig2' is deprecated; do specify 'sigma' (= sqrt(sig2))
## instead

## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
## Warning: Zero sample variances detected, have been offset away from zero

## Warning: The use of 'sig2' is deprecated; do specify 'sigma' (= sqrt(sig2))
## instead

## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing

## Warning: The use of 'sig2' is deprecated; do specify 'sigma' (= sqrt(sig2))
## instead

## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing

## Warning: The use of 'sig2' is deprecated; do specify 'sigma' (= sqrt(sig2))
## instead

## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
## Warning: Zero sample variances detected, have been offset away from zero

## Warning: The use of 'sig2' is deprecated; do specify 'sigma' (= sqrt(sig2))
## instead

## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing

## Warning: The use of 'sig2' is deprecated; do specify 'sigma' (= sqrt(sig2))
## instead

## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
## Warning: Zero sample variances detected, have been offset away from zero

## Warning: The use of 'sig2' is deprecated; do specify 'sigma' (= sqrt(sig2))
## instead

## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
## Warning: Zero sample variances detected, have been offset away from zero

## Warning: The use of 'sig2' is deprecated; do specify 'sigma' (= sqrt(sig2))
## instead

## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing

## Warning: The use of 'sig2' is deprecated; do specify 'sigma' (= sqrt(sig2))
## instead

## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
## Warning: Zero sample variances detected, have been offset away from zero

## Warning: The use of 'sig2' is deprecated; do specify 'sigma' (= sqrt(sig2))
## instead

## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing

## Warning: The use of 'sig2' is deprecated; do specify 'sigma' (= sqrt(sig2))
## instead

## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
## Warning: Zero sample variances detected, have been offset away from zero

## Warning: The use of 'sig2' is deprecated; do specify 'sigma' (= sqrt(sig2))
## instead

## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing

## Warning: The use of 'sig2' is deprecated; do specify 'sigma' (= sqrt(sig2))
## instead

## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing

## Warning: The use of 'sig2' is deprecated; do specify 'sigma' (= sqrt(sig2))
## instead

## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
## Warning: Zero sample variances detected, have been offset away from zero

## Warning: The use of 'sig2' is deprecated; do specify 'sigma' (= sqrt(sig2))
## instead

## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
## Warning: Zero sample variances detected, have been offset away from zero

## Warning: The use of 'sig2' is deprecated; do specify 'sigma' (= sqrt(sig2))
## instead

## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
## Warning: Zero sample variances detected, have been offset away from zero

## Warning: The use of 'sig2' is deprecated; do specify 'sigma' (= sqrt(sig2))
## instead

## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing

## Warning: The use of 'sig2' is deprecated; do specify 'sigma' (= sqrt(sig2))
## instead

## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
## Warning: Zero sample variances detected, have been offset away from zero

## Warning: The use of 'sig2' is deprecated; do specify 'sigma' (= sqrt(sig2))
## instead

## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
## Warning: Zero sample variances detected, have been offset away from zero

## Warning: The use of 'sig2' is deprecated; do specify 'sigma' (= sqrt(sig2))
## instead

## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing

## Warning: The use of 'sig2' is deprecated; do specify 'sigma' (= sqrt(sig2))
## instead

## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing

## Warning: The use of 'sig2' is deprecated; do specify 'sigma' (= sqrt(sig2))
## instead

## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
## Warning: Zero sample variances detected, have been offset away from zero

## Warning: The use of 'sig2' is deprecated; do specify 'sigma' (= sqrt(sig2))
## instead

## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing

## Warning: The use of 'sig2' is deprecated; do specify 'sigma' (= sqrt(sig2))
## instead

## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
## Warning: Zero sample variances detected, have been offset away from zero

## Warning: The use of 'sig2' is deprecated; do specify 'sigma' (= sqrt(sig2))
## instead

## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
## Warning: Zero sample variances detected, have been offset away from zero

## Warning: The use of 'sig2' is deprecated; do specify 'sigma' (= sqrt(sig2))
## instead

## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing

## Warning: The use of 'sig2' is deprecated; do specify 'sigma' (= sqrt(sig2))
## instead

## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
## Warning: Zero sample variances detected, have been offset away from zero

## Warning: The use of 'sig2' is deprecated; do specify 'sigma' (= sqrt(sig2))
## instead

## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing

## Warning: The use of 'sig2' is deprecated; do specify 'sigma' (= sqrt(sig2))
## instead

## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
## Warning: Zero sample variances detected, have been offset away from zero

## Warning: The use of 'sig2' is deprecated; do specify 'sigma' (= sqrt(sig2))
## instead

## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing

## Warning: The use of 'sig2' is deprecated; do specify 'sigma' (= sqrt(sig2))
## instead

## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing

## Warning: The use of 'sig2' is deprecated; do specify 'sigma' (= sqrt(sig2))
## instead

## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing

## Warning: The use of 'sig2' is deprecated; do specify 'sigma' (= sqrt(sig2))
## instead

## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing

## Warning: The use of 'sig2' is deprecated; do specify 'sigma' (= sqrt(sig2))
## instead

## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
## Warning: Zero sample variances detected, have been offset away from zero

## Warning: The use of 'sig2' is deprecated; do specify 'sigma' (= sqrt(sig2))
## instead

## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing

## Warning: The use of 'sig2' is deprecated; do specify 'sigma' (= sqrt(sig2))
## instead

## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing

## Warning: The use of 'sig2' is deprecated; do specify 'sigma' (= sqrt(sig2))
## instead

## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing

## Warning: The use of 'sig2' is deprecated; do specify 'sigma' (= sqrt(sig2))
## instead

## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing

## Warning: The use of 'sig2' is deprecated; do specify 'sigma' (= sqrt(sig2))
## instead

## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing

## Warning: The use of 'sig2' is deprecated; do specify 'sigma' (= sqrt(sig2))
## instead

## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing

## Warning: The use of 'sig2' is deprecated; do specify 'sigma' (= sqrt(sig2))
## instead

## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
## Warning: Zero sample variances detected, have been offset away from zero

## Warning: The use of 'sig2' is deprecated; do specify 'sigma' (= sqrt(sig2))
## instead

## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing

## Warning: The use of 'sig2' is deprecated; do specify 'sigma' (= sqrt(sig2))
## instead

## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing

## Warning: The use of 'sig2' is deprecated; do specify 'sigma' (= sqrt(sig2))
## instead

## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing

## Warning: The use of 'sig2' is deprecated; do specify 'sigma' (= sqrt(sig2))
## instead

## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
## Warning: Zero sample variances detected, have been offset away from zero

## Warning: The use of 'sig2' is deprecated; do specify 'sigma' (= sqrt(sig2))
## instead

## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing